PD-1 Responsive Glioblastoma Case Study: Differential Gene Expression & Pathway Analysis

Note: Some of the figures in this document use grid.arrange or ggarrange. The heights/widths arguments in grid.arrange() and ggarrange() may have to be tweaked to get the proportions of the panels to look less lopsided.

1 Setup

2 DEG analysis

2.1 Create design matrix

Create design matrix using all permunations of sample codes. 1 and 0 corresponds to which side of the comparison the sample in the sample_names column is. The sample names will just correspond to the rownames of the phenotype table. this gets combined in a for-loop later.

This design matrix is referred to as the “pooled comparison strategy” used in the analysis, as it allows us to look at the differences between different regions without resorting entirely to 1v1 comparisons. This way, we can identify gene signatures unique to a particular region compared to the rest, or to identify gene signatures that may be associated with samples containing similar clonal architectures from the previous clonevol analysis.

d1 <- ((combinat::permn(as.character(c('1', '1', '0', '0')))))
d1 <- t( unique( data.frame( matrix(unlist(d1), nrow = factorial(4), byrow = T))))
class(d1) <- "numeric"

d2 <- ((combinat::permn(as.character(c('1', '1', '1', '0')))))
d2 <- t( unique( data.frame( matrix(unlist(d2), nrow = factorial(4), byrow = T))))
class(d2) <- "numeric"

design <- cbind(d1, d2)
colnames(design) <- NULL
rownames(design) <- NULL

# clean up intermediaries
rm(d1, d2)

2.2 Generate DEG results using pooled comparisons in DESeq2

comparison_list <- NULL # store the comparisons here.
deg_list <- NULL # store the results of the comparison here.
# loop through the combinations in the design matrix to do permuted DEG tests in pooled
# comparisons using DESeq2.
for (i in 1:ncol(design)){

    sample_code <- design[, i]
    ph <- data.frame(sample_name=rownames(case.phenotype), sample_code=sample_code)

    dds <- DESeqDataSetFromMatrix(countData = gcounts, colData = ph, design = ~sample_code)
    dds <- DESeq(dds)
    res <- results(dds)
    resOrdered <- as.data.frame(res[order(res$padj),])
    resOrdered$ensid <- rownames(resOrdered)
    deg_list[[i]] <- resOrdered
    comparison_list[[i]] <- ph
}

# clean up intermediaries
rm(resOrdered, dds, res, ph, sample_code)

2.3 Clean up sample labels for visualization

# now take the comparison list into something easier to coerce into a sample v sample string
# for the volcano plot labels
comparison_names <- NULL

for (i in 1:length(comparison_list)){
    comparison <- comparison_list[[i]]
    s0 <- paste(as.character(comparison[comparison$sample_code == 0,]$sample_name), collapse = ', ')
    s1 <- paste(as.character(comparison[comparison$sample_code == 1,]$sample_name), collapse=', ')

    deg_list[[i]]$sample_comparison <- as.character(paste(s0, s1, sep=' - vs - '))
    
    comparison_names[i] <- as.character(paste(s0, s1, sep=' - vs - '))
}

names(deg_list) <- comparison_names
rm(s0, s1)

# get top 25 p values and top/bottom logFC values
top_genes <- lapply(deg_list, function(x){
    x <- x %>%  filter(!is.na(padj)) %>%
        filter(padj < 0.5) %>%
        left_join(genes, by=c('ensid'))
    x1 <- x %>% top_n(15, wt=as.numeric(-log(padj)))
    x2 <- x %>% filter(log2FoldChange > 1) %>% top_n(10, wt=abs(log2FoldChange))
    x3 <- x %>% filter(log2FoldChange < -1) %>% top_n(10, wt=abs(log2FoldChange))
    x <- unique(rbind(x1, x2, x3))

    rm(x1, x2, x3)
    return(x)
})

top_genes <- do.call(rbind, top_genes)

2.4 Remove mirrored comparisons

# coerce results to data frame and remove mirrored comparisons
mirrored_comparisons <- c("C, P - vs - A, B", "B, C - vs - A, P", "B, P - vs - A, C")
deg_df <- do.call(rbind, deg_list)
deg_df <- deg_df %>% as_tibble() %>%
    filter(!is.na(padj)) %>%
    filter(!sample_comparison %in% mirrored_comparisons) %>%
    left_join(genes, by=c('ensid'))

top_genes <- top_genes %>%
  filter(!is.na(padj)) %>%
  filter(!sample_comparison %in% mirrored_comparisons)

summary(as.factor(top_genes$sample_comparison))

## A - vs - B, C, P A, B - vs - C, P A, C - vs - B, P A, P - vs - B, C 
##                1               33               17                1 
## B - vs - A, C, P C - vs - A, B, P P - vs - A, B, C 
##               35               32               34

2.5 Visualize with volcano plot

volcano_plot_filt <- ggplot(data=deg_df[! deg_df$sample_comparison %in% c('A - vs - B, C, P', 'A, P - vs - B, C'),]) + aes(x=log2FoldChange, y=-log(padj), color=biotype) +
    geom_point() +
    geom_label_repel(data=top_genes[! top_genes$sample_comparison %in% c('A - vs - B, C, P', 'A, P - vs - B, C'),], aes(label=gene), show.legend=FALSE,  force=3, size=4) +
    facet_wrap(~ sample_comparison, scales='free',nrow=1) +
    theme_classic() + theme

volcano_plot_filt

2.6 Get GSEA scores for MSigDB oncogenic signatures

First coerce the deg results from DEseq2 into a list of ranks for each comparison using the method described here: https://stephenturner.github.io/deseq-to-fgsea/

Then, we will use fgsea to perform gene set enrichment on all the MSigDB oncogenic signatures.

deg_ranks <- lapply(deg_list, function(x){
    x %>%
        left_join(genes, by=c('ensid')) %>%
        na.omit() %>%
        distinct() %>%
        group_by(gene) %>%
        summarise(stat=mean(stat)) %>%
        deframe()
    } )

deg_ranks <- deg_ranks[! names(deg_ranks) %in% mirrored_comparisons]
deg_ranks <- deg_ranks[! names(deg_ranks) %in% c("A - vs - B, C, P","A, P - vs - B, C")]

load('../data/msigdb_gene_signatures.rdata')

Generate the enrichment scores & coerce to dataframe

# compute gene set enrichments
all_fgsea <- list(hallmark_fgsea = gsea_lister(rank_list = deg_ranks, sig_list = hallmark_gsigs, sig_name='HALLMARK'),
                  kegg_fgsea = gsea_lister(rank_list = deg_ranks, sig_list = kegg_gsigs, sig_name = 'KEGG'),
                  reactome_fgsea = gsea_lister(rank_list = deg_ranks, sig_list = reactome_gsigs, sig_name = 'REACTOME'),
                  biocarta_fgsea = gsea_lister(rank_list = deg_ranks, sig_list = biocarta_gsigs, sig_name = 'BIOCARTA')
)

Select the top pathways, filter irrelevant terms, and order the pathways via hierarchical clustering.

# select top pathways & order pathways via hclust
data_list <- lapply(all_fgsea, function(x){
  sig_name <- levels(as.factor(x$signature_name))
  x$pathway <- gsub(paste0(sig_name, '_'), '', x=x$pathway)
  x$sample_comparison <- factor(x$sample_comparison, levels=names(deg_ranks))
  
  df <- x %>%
  filter(padj < 0.05) %>% arrange(padj) %>% group_by(sample_comparison) %>%
  top_n(10, wt=-padj) %>% top_n(10, wt=abs(NES)) %>%
  arrange(abs(NES), -padj) 

  tmp <- df %>% dplyr::select(sample_comparison, pathway, NES) %>%
  group_by(sample_comparison) %>% tidyr::spread(pathway, NES) %>% 
  as.data.frame()
  
  names <- tmp$sample_comparison
  tmp$sample_comparison <- NULL
  tmp <- as.data.frame(lapply(tmp, function(x){
    x[is.na(x)] <- 0
    return(x) }))
  
  rownames(tmp) <- names
  order <- hclust( dist(t(tmp), method = "euclidean"), method = "ward.D" )$order
  df$pathway <- factor(df$pathway,levels=colnames(tmp)[order])


  df <- df[!is.na(df$pathway),]
  df <- df[df$pathway != 'NA',]

  rm(names, tmp)

  return(df)
})

terms <- c('INFLUENZA', 'CELL_CYCLE', 'MITOTIC', 'S_PHASE', 'G1_S_TRANSITION', 'ASTHMA', 'AUTOIMMUNE_THYROID_DISEASE',
'INFECTION','DISEASE', 'INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION','METABOLISM', 'APICAL_SURFACE','GENESIS',
 '43S', 'SYNAP','NEUROTRANSMITTER', 'GABA', 'GLUTAMATE', 'POTASSIUM_CHANNELS', 'CHANNEL_TRANSPORT', 'CELL_LINEAGE',
 'LUPUS', 'NEURONAL',  'RP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE', 'TRANSLATION', 'proteasome', 'complement', 'homologous_recombination',
"peptide_chain",'second_messenger', 'REACTIVE_OXYGEN_SPECIES', 'NONSENSE_MEDIATED_DECAY', 'G2M_checkpoint', 'G2_M_CHECKPOINT', 'IMMUNOREGULATORY_INTERACTIONS_BETWEEN',
'RIBOSOME', 'TIGHT_JUNCTION' ) 


data_filt <- lapply(data_list, function(x){

  y <- x
  for (term in terms){
    y <- y %>% filter(!grepl(term, x=pathway, ignore.case=TRUE))
  }

df <- y

tmp <- df %>% dplyr::select(sample_comparison, pathway, NES) %>%
  group_by(sample_comparison) %>% tidyr::spread(pathway, NES) %>% 
  as.data.frame()
  
  names <- tmp$sample_comparison
  tmp$sample_comparison <- NULL
  tmp <- as.data.frame(lapply(tmp, function(x){
    x[is.na(x)] <- 0
    return(x) }))
  
  rownames(tmp) <- names
  order <- hclust( dist(t(tmp), method = "euclidean"), method = "ward.D" )$order
  df$pathway <- factor(df$pathway,levels=colnames(tmp)[order])


  df <- df[!is.na(df$pathway),]
  df <- df[df$pathway != 'NA',]

  rm(names, tmp)

  return(df)

})

2.7 Visualize gene set enrichments as a heatmap

Make one plot per gene set database and visualize as a set of stacked heatmaps.

# set up themes for ggplotting
theme <- theme(text = element_text(size=16),
               axis.text.y = element_text(size=18, face='bold'),
               axis.text.x = element_text(size=18, face='bold'),
               axis.title = element_text(size=24, face = "bold"),
               legend.text= element_text(size=18),
               legend.title= element_text(size=24, face='bold'),
               strip.text=element_text(size=24, face='bold'),
               plot.title=element_text(size=28, face='bold'))

plot_list <- lapply(data_filt, function(x){
  plot <- ggplot(data=x) +
  aes(x=sample_comparison, y=pathway, fill=as.numeric(NES)) +
  geom_tile(color='white', size=1) +
  theme_classic() +
  theme + theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  scale_x_discrete(drop=FALSE) +
  labs(x ='', y=levels(as.factor(x$signature_name))) +
  scale_fill_gradient2(low = "lightseagreen", mid='white', high = "tomato", name="NES")  +
  guides(fill = guide_colourbar(barwidth = 1.5, barheight = 25)) 
  return(plot)
})

plot_list[[4]] <- plot_list[[4]] + theme_classic() + theme + rotate_x_text(angle=45)

names(plot_list) <- c('HALLMARK', 'KEGG', 'REACTOME', 'BIOCARTA')
pathway_hm <- ggarrange(plotlist=plot_list, nrow=length(all_fgsea), align='v',  common.legend=TRUE, legend='right', heights=c(1.2,0.6,0.5,1.6))

pathway_hm

3 Integrate TCGA expression information

3.1 PCA

case_int_pcs <- ind

surv_pca <- ggplot(data = case_int_pcs) + 
  aes(x=Dim.1, y=Dim.2, color=log(MonthsOverallSurvival), label=label, shape=VitalStatusFactor) +
  geom_density2d(color='black',alpha=0.5, bins=50, aes(shape=NULL)) + 
  geom_point(alpha=0.8, size=6) +
  geom_label_repel(size=8, show.legend=FALSE, force=10, point.padding = unit(0.5, 'lines'), fontface='bold') +
  labs(x='PC-1: 14.3% of variance explained', y='PC-2: 9% of variance explained') +
  theme_classic() + 
  theme(plot.margin=margin(0.5,0.5,0.5,0.5, "in")) +
  theme + scale_color_gradient(low='ivory', high='black', limits=c(-2,5), breaks=seq(-2,5,length.out = 5)) + 
  guides(color = guide_colourbar(barwidth = 1.5, barheight = 15)) 

surv_pca

## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 151 rows containing missing values (geom_label_repel).

3.2 With distributions of specific checkpoint and immune genes

Objective: get a measure (a p-value) of whether the variance between my four samples is greater than if you were to take the variance four random TCGA samples. To do this, do resamplings of 10% of the total possible combinations of 4 tcga-GBM samples from a pool of 155 samples.

checkpoint.genes <- c("PD1","PDL1","PDL2","TIM3","CD28","CD3","CD4","FOXP3","DCLAMP")
val_cols <- checkpoint.genes
chk <- ind
chk$submitter_id <- as.character(chk$submitter_id)

# separate out samples from case study
ref_df <- chk[1:4,]
ref_df <- data.frame(var=sapply(ref_df[, val_cols], var),
                      samp_type=as.character('test'),
                      val=val_cols)

set.seed(0)

length(combn(chk[- c(1:4),]$submitter_id, 4, simplify=FALSE))
boot_pool <- chk[- c(1:4),]
# grab all unique combinations of 4 from the pool of 155 TCGA samples
# and then only take 10% of that cause there's too many
combinations <- combn(boot_pool$submitter_id, 4, simplify=FALSE)
combinations <- base::sample(combinations, length(combinations)*0.1)

# function to take variance of each combination of 4 for each gene/column
boot_fx <- function(x){
  # grab the combination of 4
  x <- boot_pool %>%
       filter(submitter_id %in% x) %>%
       dplyr::select(val_cols)
  # take the variance for each gene 
  x <- data.frame(var=sapply(x, var),
                  val=val_cols,
                  samp_type=as.character('random_boot'))
  return(x)
}

# do the resampling; using 3 cores (but the more the merrier!), and then bind the rows.
boot_pool <- bind_rows(pbmclapply(combinations, boot_fx, mc.cores=2))
all_boots <- bind_rows(ref_df, boot_pool)

Get a measure of whether the variance between my four samples is greater than if you were to take four random TCGA samples.

pvals <- lapply(levels(factor(all_boots$val)), function(i){
    x <- all_boots %>% filter(val == i) %>% as.data.frame()
    test_point <- x %>% filter(samp_type == "test") %>% dplyr::select(var) %>% deframe()
    z <- stats::density(x$var, n=length(x$var))
    xval <- as.numeric(as.matrix(as.numeric(z$x)))
    yval <- as.numeric(as.matrix(as.numeric(z$y)))

    # test right-sided
    rpval <- auc(x=xval, y =yval, from=test_point, to=max(xval) )
    # test left-sided
    lpval <- auc(x=xval, y =yval, from=min(xval), to=test_point )
    # pick the smallest of the two left-sided and right-sided tests
    # to get the auc (p-value) for the closest end of the distribution
    result <- min(c(rpval, lpval))
    return(result)
})

names(pvals) <- levels(factor(all_boots$val))
pvals <- unlist(pvals)

pvaldf <- data.frame(pval=pvals,
                    gene=names(pvals))

3.3 plot the checkpoint violins

ind2 <- tidyr::gather(chk, gene, cpm, val_cols)
ind2$gene <- as.factor(ind2$gene)
ind2$cpm <- as.numeric(ind2$cpm)


ind2 <- ind2 %>% inner_join(pvaldf, by=c('gene'))
ind2$sample <- as.factor(ind2$sample)

Note: we saved a pre-computed version of this to reduce runtime and will use it in this RMarkdown document.

checkpoint_violins <- ggplot(data=ind2) + aes(y=cpm, x=gene, color=sample) + 
  #geom_violin(alpha=0.05, aes(color=NULL), fill="#7f7f7f", show.legend = FALSE) + 
  #geom_boxplot(width=0.25, aes(color=NULL, fill=NULL), show.legend = FALSE ) +
  geom_beeswarm(data=ind2[ind2$sample =='TCGA-GBM',], alpha=0.4, size=2, show.legend=FALSE) +
  geom_point(data=ind2[ind2$sample !='TCGA-GBM',], size=6, alpha=0.75) +
  theme_classic() + theme + labs(y='Normalized Expression', x='', color='Sample') + 
  # don't add duplicate labels
  geom_text(data=ind2[(ind2$pval<0.05) & (ind2$sample=='Primary'),],
            y=-5,size=6, show.legend=FALSE,
            aes(x=gene, label=paste('P =',round(pval, digits=4)), color=NULL)) +
            scale_color_manual(values = cbpal, aesthetics = 'color')

checkpoint_violins

Supplementary Figure 3 - overlay onto the PCA

gene_pca <- ggplot(data = ind2) + 
  aes(x=Dim.1, y=Dim.2, color=cpm, label=label) +
  geom_density2d(color="#7f7f7f", bins=50, alpha=0.7) + 
  geom_point(alpha=0.7, size=5) +
  geom_label_repel(size=8, show.legend=FALSE, force=10, min.segment.length = 10, segment.size = 1) +
  labs(x='PC-1: 14.3% of variance explained', y='PC-2: 9% of variance explained') +
  facet_wrap(~ gene) +
  theme_classic() +
  theme + scale_color_gradient(low="ivory", high="black", name='Norm. Exprs.', limits=c(-5,10), breaks=seq(-5, 10, by=5)) +
  guides(color = guide_colourbar(barwidth = 1.5, barheight = 25)) 

# pdf('../../figures/tcga_ext.pdf', height=20, width=24)
gene_pca

# dev.off()

4 Manuscript Figure 4

tcga_integration2 <- ggarrange(surv_pca, checkpoint_violins, nrow=2, heights=c(2,0.7), align='v')

## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 151 rows containing missing values (geom_label_repel).

#pdf('../../figures/fig4.pdf', height=18, width=34)
grid.arrange(pathway_hm, tcga_integration2, ncol=2, widths=c(0.7, 1.5))

#dev.off()

5 Survival integration

surv_dens <- ggplot(data = case_int_pcs) + 
    aes(x=MonthsOverallSurvival) + 
    geom_density(fill='gray', alpha=0.7, show.legend = FALSE) + 
    theme_classic() +
    theme + 
    geom_vline(data=case_int_pcs[!is.na(case_int_pcs$label),], size=1.5, color='blue',
               aes(xintercept=MonthsOverallSurvival), show.legend = FALSE, linetype='dashed') + 
    geom_text(data=case_int_pcs[1,], 
              aes(x = MonthsOverallSurvival, y=0.02), label = 'Patient Survival', color='blue',
                    show.legend = FALSE,angle=90, size=8, nudge_x =3, fontface='bold') + 
    scale_y_continuous(limits = c(0,0.05)) 

#pdf('../../figures/survival_density_tcga.pdf', height=4, width=6)
surv_dens

## Warning: Removed 4 rows containing non-finite values (stat_density).

#dev.off()

The patient survived almost 1 and a half times longer than the third quartile of tcga survival.

summary(case_int_pcs$MonthsOverallSurvival)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.16    5.08   11.84   14.15   17.75   88.14       4

Use ECDF on the TCGA cohort to find what the probability is that our patient survived this long by chance.

ecdf_surv <- ecdf(case_int_pcs$MonthsOverallSurvival)

plot(ecdf_surv)

paste0("The probability that the patient would have survived 25 months by chance alone is ", round((1 - ecdf_surv(25))*100, 2), '%' )

## [1] "The probability that the patient would have survived 25 months by chance alone is 15.89%"

6 GSEA for Zhao signatures

Taking the pre-treatment and post-treatment responder and non-responder signatures in figure 3b, and ext data fig 7, do single sample GSEA on the relevant gene sets to see if there is any association.

set.seed(0)

zhao_enrichments <- list(pretx_responder_down = gsea_lister(rank_list = deg_ranks, sig_list =zhao_responder_pretx_dn,
                                                            sig_name='PRE-TX_DN'),
                  pretx_responder_up= gsea_lister(rank_list = deg_ranks, sig_list = zhao_responder_pretx_up, sig_name = 'PRE-TX_UP'),
                  posttx_responder_down = gsea_lister(rank_list = deg_ranks, sig_list = zhao_responder_posttx_dn, sig_name = 'POST-TX_DN'),
                  posttx_responder_up = gsea_lister(rank_list = deg_ranks, sig_list = zhao_responder_posttx_dn, sig_name = 'POST-TX_UP')
)

zhao_enrichments %>% bind_rows() %>% filter(padj < 0.05)

## # A tibble: 26 x 7
##    pathway               pval   padj   NES  size sample_comparis… signature_name
##    <chr>                <dbl>  <dbl> <dbl> <int> <chr>            <chr>         
##  1 XU_GH1_EXOGENOUS_… 0.00149 0.0104  1.80    92 A, B - vs - C, P PRE-TX_DN     
##  2 GSE21678_WT_VS_FO… 0.00140 0.0104  1.54   162 A, B - vs - C, P PRE-TX_DN     
##  3 GSE17721_LPS_VS_P… 0.00121 0.0169  1.62   180 P - vs - A, B, C PRE-TX_DN     
##  4 FAELT_B_CLL_WITH_… 0.0168  0.0454 -1.57    48 C - vs - A, B, P PRE-TX_DN     
##  5 MARCHINI_TRABECTE… 0.0235  0.0454 -1.67    18 C - vs - A, B, P PRE-TX_DN     
##  6 GSE17721_LPS_VS_P… 0.0264  0.0454 -1.37   184 C - vs - A, B, P PRE-TX_DN     
##  7 XU_GH1_EXOGENOUS_… 0.0163  0.0454 -1.52    91 C - vs - A, B, P PRE-TX_DN     
##  8 MTOR_UP.N4.V1_DN   0.0213  0.0454 -1.37   155 C - vs - A, B, P PRE-TX_DN     
##  9 GSE45739_NRAS_KO_… 0.00256 0.0359 -1.48   174 C - vs - A, B, P PRE-TX_DN     
## 10 GSE21678_WT_VS_FO… 0.0292  0.0454 -1.34   160 C - vs - A, B, P PRE-TX_DN     
## # … with 16 more rows

znmed_list <- lapply(zhao_enrichments, function(x){
  x$sample_comparison <- factor(x$sample_comparison, levels=names(deg_ranks))
  df <- x %>% filter(padj < 0.05)
  
  plot <- ggplot(data=df) +
  aes(x=sample_comparison, y=pathway, fill=as.numeric(NES)) +
  geom_tile(color='white', size=1) +
  theme_classic() +
  theme + theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  scale_x_discrete(drop=FALSE) +
  labs(x ='', y=levels(as.factor(x$signature_name))) +
  scale_fill_gradient2(low = "lightseagreen", mid='white', high = "tomato", name="NES")  +
  guides(fill = guide_colourbar(barwidth = 1.5, barheight = 25)) 
  return(plot)
})

znmed_list[[4]] <- znmed_list[[4]] + theme_classic() + theme + rotate_x_text(angle=45)

zhao_hm <- ggarrange(plotlist=znmed_list, nrow=length(znmed_list), align='v',  common.legend=TRUE, legend='right', heights=c(4,1,1,1.7))

zhao_hm

#pdf('../../figures/Sfig3_tcga_zhao_ext.pdf', height=18, width=34)
grid.arrange(zhao_hm, gene_pca, ncol=2, widths=c(1, 1.5))

## Warning: Removed 1359 rows containing missing values (geom_label_repel).

#dev.off()

7 RSession Info

sessionInfo()

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] lattice_0.20-38             ggbeeswarm_0.7.0           
##  [3] fgsea_1.10.1                Rcpp_1.0.3                 
##  [5] ggpubr_0.2.4                magrittr_1.5               
##  [7] org.Hs.eg.db_3.8.2          AnnotationDbi_1.46.1       
##  [9] gridExtra_2.3               ggrepel_0.8.1              
## [11] combinat_0.0-8              forcats_0.4.0              
## [13] stringr_1.4.0               dplyr_0.8.3                
## [15] purrr_0.3.3                 readr_1.3.1                
## [17] tidyr_1.0.0                 tibble_2.1.3               
## [19] ggplot2_3.2.1               tidyverse_1.2.1            
## [21] edgeR_3.26.8                limma_3.40.6               
## [23] DESeq2_1.24.0               SummarizedExperiment_1.14.1
## [25] DelayedArray_0.10.0         BiocParallel_1.18.1        
## [27] matrixStats_0.55.0          Biobase_2.44.0             
## [29] GenomicRanges_1.36.1        GenomeInfoDb_1.20.0        
## [31] IRanges_2.18.3              S4Vectors_0.22.1           
## [33] BiocGenerics_0.30.0        
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       ggsignif_0.6.0         htmlTable_1.13.2      
##  [4] XVector_0.24.0         base64enc_0.1-3        rstudioapi_0.10       
##  [7] bit64_0.9-7            fansi_0.4.0            lubridate_1.7.4       
## [10] xml2_1.2.2             splines_3.6.1          geneplotter_1.62.0    
## [13] knitr_1.26             zeallot_0.1.0          Formula_1.2-3         
## [16] jsonlite_1.6           broom_0.5.2            annotate_1.62.0       
## [19] cluster_2.1.0          compiler_3.6.1         httr_1.4.1            
## [22] backports_1.1.5        assertthat_0.2.1       Matrix_1.2-17         
## [25] lazyeval_0.2.2         cli_1.1.0              acepack_1.4.1         
## [28] htmltools_0.4.0        tools_3.6.1            gtable_0.3.0          
## [31] glue_1.3.1             GenomeInfoDbData_1.2.1 fastmatch_1.1-0       
## [34] cellranger_1.1.0       vctrs_0.2.0            nlme_3.1-142          
## [37] xfun_0.11              rvest_0.3.5            lifecycle_0.1.0       
## [40] XML_3.98-1.20          MASS_7.3-51.4          zlibbioc_1.30.0       
## [43] scales_1.0.0           hms_0.5.2              RColorBrewer_1.1-2    
## [46] yaml_2.2.0             memoise_1.1.0          rpart_4.1-15          
## [49] latticeExtra_0.6-28    stringi_1.4.3          RSQLite_2.1.2         
## [52] genefilter_1.66.0      checkmate_1.9.4        rlang_0.4.1           
## [55] pkgconfig_2.0.3        bitops_1.0-6           evaluate_0.14         
## [58] labeling_0.3           htmlwidgets_1.5.1      cowplot_1.0.0         
## [61] bit_1.1-14             tidyselect_0.2.5       R6_2.4.1              
## [64] generics_0.0.2         Hmisc_4.3-0            DBI_1.0.0             
## [67] pillar_1.4.2           haven_2.2.0            foreign_0.8-72        
## [70] withr_2.1.2            survival_3.1-7         RCurl_1.95-4.12       
## [73] nnet_7.3-12            modelr_0.1.5           crayon_1.3.4          
## [76] utf8_1.1.4             rmarkdown_1.17         locfit_1.5-9.1        
## [79] readxl_1.3.1           data.table_1.12.6      blob_1.2.0            
## [82] digest_0.6.22          xtable_1.8-4           munsell_0.5.0         
## [85] beeswarm_0.2.3         vipor_0.4.5